# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# REPLICATION MATERIAL FOR:
# Who is afraid of the Istanbul Convention? 
# Explaining opposition to and support for gender equality
# 
# Michael C. Zeller
# with Andrea Krizsan and Conny Roggeband
#
# 2024
#
# R version 4.3.1
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

setwd("~/Downloads/Supplemental material")

library(dplyr)
library(ggplot2)
library(reshape2)
library(QCA)
library(SetMethods)
library(readxl)
library(ComplexUpset)
library(eulerr)
if (!require("xtable")) install.packages("xtable")
library(kableExtra)
library(knitr)

ICdata <- read.csv("database_IC-qca.csv", row.names=1, header = TRUE)

# Outcome -----------------------------------------------------------------

# IC (fuzzy) = full Istanbul Convention ratification
#### 1 - signed and ratified; 
#### 0.75 - signed and ratified with declaration; 
#### 0.25 - signed but not ratified; 
#### 0 - not signed/denounced

# Conditions --------------------------------------------------------------

# LDB - liberal democratic backsliding, using v2xlibdem variable (2021 score subtracted by 2013 score)
#### Calibration (direct method)
#### Robustness test w/ EDB, egalitarian democratic backsliding, using v2xegaldem variable (2021 score subtracted by 2013 score)
#### Robustness test w/ SLD, strong liberal democracy, using v2xlibdem variable (2021 score)
#### Robustness test w/ SED, strong egalitarian democracy, using v2xegaldem variable (2021 score)

# RGOV (fuzzy) = solid right-wing government (through signature and ratification phase) 
#### Calibration (direct assignment): 
#### 1    - solid right throughout signature and ratification; 
#### 0.75 - not solid right at signature, solid right since signature; 
#### 0.45 - solid right at signature, not solid right in post-signature; 
#### 0.25 - not solid right signature - left right changes afterwards; 
#### 0    - not solid right, neither at signature nor ratification

# OPPG - opposition to 'gender ideology'
#### Calibration (direct assignment): derived from unweighted index of EVS/ESS data,
#### A124_09 (EVS: homosexual neighbours are undesirable), 
#### D081_N (EVS: disagree that homosexuals are good parents), 
#### F118_N (EVS: homosexuality is never justifiable))

# HRM (crisp) = homogenous religious majority
#### if a country is predominantly Catholic, Orthodox, or Islamic (EVS: F025_MAJ), 
#### and if that predominance is over 50% of the population (A006), 
#### then the country is scored as 1; otherwise scored as 0

# WPOW - women's political empowerment: derived from unweighted index of V-Dem data
#### v2x_genpp_2020 (women’s political empowerment women's political participation) and
#### v2x_cspart (civil society participation)

# WEQU - strong social support for women’s equality
#### Calibration (direct assignment): derived from EVS/ESS data
#### C001 (EVS: men should have more right to a job than women)

# TIDY THE DATA

# removing rows
row.names.remove <- c("Andorra", "Liechtenstein", "Luxembourg", "Malta", "Monaco", "Moldova", "San Marino")

ICdata <- ICdata[!(row.names(ICdata) %in% row.names.remove), ]

ICdatax <- ICdata[!(row.names(ICdata) %in% row.names.remove), ]

# removing excess columns
ICdatax <- ICdatax[-c(2:7,24:106)]

# i. calibration ----------------------------------------------------------

# LDB: liberal democratic backsliding

# ICdatax$LDEMBACK <- ICdatax$v2x_libdem_2020 - ICdatax$v2x_libdem_2013
ICdatax$LDEMBACK <- ICdatax$v2x_libdem_2021 - ICdatax$v2x_libdem_2013

ICdatax$LDB <- calibrate(ICdatax$LDEMBACK,
                         type = "fuzzy",
                         thresholds = c(0, -0.08, -0.1))
ICdatax$LDB

# WPOW: women's political empowerment 
ICdatax$WOMPOW <- (ICdatax$v2x_gender_2020 + ICdatax$v2x_cspart_2020)/2

ICdatax$WPOW <- calibrate(ICdatax$WOMPOW,
                          type = "fuzzy",
                          thresholds = c(0.35, 0.9, 0.97))
ICdatax$WPOW

####
# alternative analysis conditions
####

# EDB: egalitarian democratic backsliding
#ICdatax$EDEMBACK <- ICdatax$v2x_egaldem_2020 - ICdatax$v2x_egaldem_2013
ICdatax$EDEMBACK <- ICdatax$v2x_egaldem_2021 - ICdatax$v2x_egaldem_2013

ICdatax$EDB <- calibrate(ICdatax$EDEMBACK,
                         type = "fuzzy",
                         thresholds = c(0.04, -0.065, -0.2))
ICdatax$EDB

# SLD: strong liberal democracy, v2xlibdem variable (2021 score)
ICdatax$SLD <- calibrate(ICdatax$v2x_libdem_2021,
                         type = "fuzzy",
                         thresholds = c(0.15, 0.775, 0.79))
ICdatax$SLD

# SED, strong egalitarian democracy, v2xegaldem variable (2021 score)
ICdatax$SED <- calibrate(ICdatax$v2x_egaldem_2021,
                         type = "fuzzy",
                         thresholds = c(0.25, 0.65, 0.85))
ICdatax$SED

# WPAR: women's political participation
ICdatax$WPAR <- calibrate(ICdatax$v2x_genpp_2020,
                          type = "fuzzy",
                          thresholds = c(0.8, 0.9945, 1))
ICdatax$WPAR



# iii. for table of all democracy conditions ------------------------------

ICdataDEM <- select(ICdatax, c("LDB","EDB","SLD","SED"))


# Right-wing govs. figure -------------------------------------------------

# Istanbul Convention adoption and solidly right-wing govs in 40 countries

govs <- read.csv("R_govs.csv")

govs$Order <- factor(govs$Order,
                       levels = c("1", "2", "3", "4", "5"))

ICgraphic <- ggplot(govs, aes(x=Country, y=Days, fill=Order)) +
  geom_bar(stat = "identity", width = 0.7, position = position_stack(reverse = T)) +
  scale_fill_manual("legend", values = c("1" = "grey50", "2" = "white", 
                                         "3" = "grey50", "4" = "white", 
                                         "5" = "grey50")) +
  scale_y_continuous("Year", breaks = seq(0,4015,365),
                     labels = paste0(c("2011","2012","2013","2014","2015",
                                       "2016","2017","2018","2019","2020",
                                       "2021","2022"))) +
  geom_point(aes(y=Signed), shape=1, size=2) +
  geom_point(aes(y=Ratified), shape=23, fill="black", size=2) +
  geom_point(aes(y=Withdrawn), shape=4, size=2) +
  coord_flip() +
  theme(
    # axis.ticks.y = element_blank(),
    # axis.text.y = element_text(hjust = 1.3),
    panel.background = element_rect(fill = NA),
    panel.ontop = TRUE,
    panel.grid.major.y = element_blank(),
    # explicitly set the horizontal lines (or they will disappear too)
    panel.grid.major.x = element_line(size=.2, color="black"),
    panel.grid.minor = element_blank(),
    legend.position = "none",
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
  ) + 
  scale_x_discrete(limits=rev)
ICgraphic + geom_vline(xintercept = ((-0.51:39.51)+1), linewidth=0.5, color = "grey80")


# Necessity ---------------------------------------------------------------

ICdatay <- select(ICdatax, c("IC","LDB","RGOV","OPPG","HRM","WPOW","WEQU",
                             "EUC"))

SUIN_I <- superSubset(data = ICdatay,
                      conditions = c("LDB","RGOV","OPPG","HRM","WPOW",
                                     "WEQU","EUC"),
                      outcome = "IC",
                      relation = "nec",
                      incl.cut = 0.9,
                      cov.cut = 0.6,
                      ron.cut = 0.5,
                      depth = 2)
# SUIN_I

SUIN_IC <- superSubset(data = ICdatay,
                       conditions = c("LDB","RGOV","OPPG","HRM","WPOW",
                                      "WEQU","EUC"),
                       outcome = "IC",
                       relation = "nec",
                       incl.cut = 0.94,
                       cov.cut = 0.6,
                       ron.cut = 0.5,
                       depth = 2)

# pimplot(data = ICdatay, results = SUIN_IC, outcome = "IC",necessity = TRUE, jitter = TRUE, font = "mono")

library(yarrr)

nec <- compute("~RGOV + ~OPPG", data = ICdatay) # compute nec. set membership scores
col <- rep(yarrr::transparent("white", trans.val = .9), nrow(ICdatay)) # define a vector of colors
col[nec > 0.5 & ICdatay$IC < 0.5] <- "grey70" # zone 4
col[nec < 0.5 & ICdatay$IC < 0.5] <- "grey70" # zone 4
col[nec < 0.5 & ICdatay$IC >= 0.5] <- "black" # zone 5
# col[nec > 0.5 & nec < 0.8 & ICdatay$IC > 0.5] <- "black"
clabels <- logical(nrow(ICdatay))
clabels[c(26)] <- TRUE # only these labels to print

XYplot("~RGOV + ~OPPG", "IC", relation = "necessity", data = ICdatay,
       jitter = TRUE, amount = 0.03, clabels = clabels,
       xlab = "~RGOV + ~OPPG (socio-political liberalism?)",
       col = "black", bg = col)



# Truth table -------------------------------------------------------------

ttIC <- truthTable(data=ICdatay, outcome = "IC",
                   conditions = c("LDB","RGOV","OPPG","HRM","WPOW","WEQU",
                                  "EUC"),
                   incl.cut=0.85,
                   sort.by=c("incl", "n"), 
                   complete=TRUE,
                   show.cases=TRUE, 
                   PRI=TRUE)
# ttIC

# stargazerTT(ttIC, show.cases = TRUE)


# Upset plot --------------------------------------------------------------

ICdataz <- subset(ICdatay, select = -c(IC))
ICdataz <- ifelse(ICdataz > 0.5, 1, 0)
ICvennZ <- as.data.frame(ICdataz)


upset(data = ICvennZ, 
      intersect = c("LDB", "RGOV", "OPPG", "HRM", "WPOW", "WEQU", "EUC"),
      base_annotations = list(
        'Intersection size'=(
          intersection_size(
            mapping=aes(fill='bars_color')) +
            scale_y_continuous(breaks=seq(0,10,1), limits=c(0,10)) +
            scale_fill_manual(values = c('bars_color'='black'), guide='none'))
      ),
      queries=list(upset_query(set='HRM', fill='black'),
                   upset_query(set='WEQU', fill='black'),
                   upset_query(set='OPPG', fill='black'),
                   upset_query(set='WPOW', fill='black'),
                   upset_query(set='LDB', fill='black'),
                   upset_query(set='EUC', fill='black'),
                   upset_query(set='RGOV', fill='black'),
                   upset_query(intersect=c('WPOW', 'WEQU'), 
                               color='black', fill='grey80', 
                               only_components=c('intersections_matrix', 'Intersection size')),
                   upset_query(intersect=c('HRM', 'OPPG', 'EUC'), 
                               color='black', fill='grey80', 
                               only_components=c('intersections_matrix', 'Intersection size')),
                   upset_query(intersect=c('LDB', 'OPPG', 'HRM'), 
                               color='black', fill='grey80', 
                               only_components=c('intersections_matrix', 'Intersection size')),
                   upset_query(intersect=c('OPPG', 'WPOW', 'WEQU'), 
                               color='black', fill='grey80', 
                               only_components=c('intersections_matrix', 'Intersection size')),
                   upset_query(intersect=c('LDB', 'HRM', 'WEQU'), 
                               color='black', fill='grey80', 
                               only_components=c('intersections_matrix', 'Intersection size')),
                   upset_query(intersect=c('LDB', 'OPPG', 'HRM', 'EUC'), 
                               color='black', fill='grey80', 
                               only_components=c('intersections_matrix', 'Intersection size')),
                   upset_query(intersect=c('LDB', 'OPPG', 'HRM', 'WPOW'), 
                               color='black', fill='grey80', 
                               only_components=c('intersections_matrix', 'Intersection size')),
                   upset_query(intersect=c('RGOV', 'OPPG', 'HRM', 'EUC'), 
                               color='black', fill='grey80', 
                               only_components=c('intersections_matrix', 'Intersection size')),
                   upset_query(intersect=c('HRM', 'WPOW', 'WEQU'), 
                               color='black', fill='grey80', 
                               only_components=c('intersections_matrix', 'Intersection size')),
                   upset_query(intersect=c('RGOV', 'WPOW', 'WEQU'), 
                               color='black', fill='grey80', 
                               only_components=c('intersections_matrix', 'Intersection size')),
                   upset_query(intersect=c('HRM', 'WPOW'), 
                               color='black', fill='grey80', 
                               only_components=c('intersections_matrix', 'Intersection size')),
                   upset_query(intersect=c('RGOV', 'HRM'), 
                               color='black', fill='grey80', 
                               only_components=c('intersections_matrix', 'Intersection size')),
                   upset_query(intersect=c('RGOV'), 
                               color='black', fill='grey50', 
                               only_components=c('intersections_matrix', 'Intersection size')),
                   upset_query(intersect=c('LDB', 'RGOV', 'OPPG'), 
                               color='black', fill='grey50', 
                               only_components=c('intersections_matrix', 'Intersection size')),
                   upset_query(intersect=c('RGOV', 'OPPG', 'WPOW'), 
                               color='black', fill='grey50', 
                               only_components=c('intersections_matrix', 'Intersection size')),
                   upset_query(intersect=c('RGOV', 'OPPG', 'HRM'), 
                               color='black', fill='grey50', 
                               only_components=c('intersections_matrix', 'Intersection size')),
                   upset_query(intersect=c('LDB', 'RGOV', 'OPPG', 'HRM', 'EUC'), 
                               color='black', fill='grey50', 
                               only_components=c('intersections_matrix', 'Intersection size'))),
      name = 'Intersections',
)


# Sufficiency solutions ---------------------------------------------------

sol_c <- minimize(ttIC, details=TRUE, show.cases=TRUE)
# sol_c

sol_i <- minimize(ttIC, details = TRUE, 
                  include = "?",
                  dir.exp = c(0, 0, 0, 0, 1, 1, 1))
# sol_i

# stargazerSol(sol_i,
#              outcome = "IC",
#              show.cases = T,
#              type = "latex",
#              title = "",
#              out = NULL,
#              digits = 3)

sol_p <- minimize(ttIC, include = "?",
                  details = TRUE, show.cases = TRUE)
# sol_p

library(yarrr)
# i_sol <- compute("~LDB*EUC + LDB*~RGOV + ~LDB*~OPPG*HRM + ~LDB*~OPPG*WPOW*WEQU +
#        ~RGOV*~HRM*WPOW*WEQU", data = ICdatay) # compute sol_i set membership scores
p_sol <- compute("WEQU + ~LDB*EUC + LDB*~RGOV + ~OPPG*HRM", data = ICdatay) # compute sol_p set membership scores
col <- rep(yarrr::transparent("white", trans.val = .9), nrow(ICdatay)) # define a vector of colors
col[p_sol < 0.5 & ICdatay$IC < 0.5] <- "black" # zone 4
col[p_sol < 0.5 & ICdatay$IC >= 0.5] <- "grey60" # zone 5
col[p_sol > 0.5 & ICdatay$IC > 0.5 & ICdatay$IC < 0.9] <- "grey60" # zone 5
# col[nec > 0.5 & nec < 0.8 & ICdatay$IC > 0.5] <- "black"
clabels <- logical(nrow(ICdatay))
clabels[c(8, 28, 9, 15, 30)] <- TRUE # only these labels to print

# TO ENSURE THE LABELS OF JITTERED POINTS ON TOP ARE VISIBLE
ICdatay$IC_ratification <- ifelse(ICdatay$IC > 0.8, ICdatay$IC-0.02, ICdatay$IC)

# TO ENSURE POLAND LABEL IS VISIBLE
# write similar ifelse function with multiple conditions?

## i_sol
# XYplot("~LDB*EUC + LDB*~RGOV + ~LDB*~OPPG*HRM + ~LDB*~OPPG*WPOW*WEQU +
#        ~RGOV*~HRM*WPOW*WEQU", IC_ratification, relation = "sufficiency", data = ICdatay, jitter = T, amount = 0.03, xlab = "Solution", clabels = clabels, cex = 0.8, col = "black", bg = col)

## p_sol
XYplot("WEQU + ~LDB*EUC + LDB*~RGOV + ~OPPG*HRM", IC_ratification, relation = "sufficiency", data = ICdatay, jitter = T, amount = 0.03, xlab = "Solution", clabels = clabels, cex = 0.8, col = "black", bg = col)


# Solution (parsimonious) terms and overlaps ------------------------------

fit <- euler(c("role \n models"=13, "EU \n dependents"=5, "pro-EU \n governments"=4, "tame \n Catholicism"=2,
               "role \n models&tame \n Catholicism"=2, "role \n models&pro-EU \n governments&tame \n Catholicism"=1))

plot(fit, fills = "transparent", 
     labels = list(font = 2, fontsize=9), 
     quantities = list(cex = 1))


# Non-outcome sufficiency solution (intermediate) -------------------------

ttNIC <- truthTable(data=ICdatay, outcome="~IC",
                    conditions = c("LDB","RGOV","OPPG","HRM","WPOW",
                                   "WEQU","EUC"),
                    incl.cut=0.75,
                    sort.by=c("incl", "n"), 
                    complete=TRUE,
                    show.cases=TRUE, 
                    PRI=TRUE)
# ttNIC

sol_Nc <- minimize(ttNIC, details=TRUE, show.cases=TRUE)
# sol_Nc

sol_Ni <- minimize(ttNIC, details = TRUE, 
                   include = "?",
                   dir.exp = c(1, 1, 1, 1, 0, 0, 0))
# sol_Ni

sol_Np <- minimize(ttNIC, include = "?", 
                   details = TRUE, show.cases = TRUE)
# sol_Np

# pimplot(data = ttNIC, results = sol_Ni, outcome = "~IC", all_labels = TRUE, jitter = TRUE, font = "mono", fontsize = 3)

# Identify contradictory simplifying assumptions: 
CSA <- LR.intersect(sol_Np, sol_p)

# Avoid contradictory simplifying assumptions: 
ttNIC_esa <- esa(oldtt = ttNIC, contrad_rows = c(CSA))

sol_Nc_esa <- minimize(ttNIC_esa, details=TRUE, show.cases=TRUE)
# sol_Nc_esa

sol_Ni_esa <- minimize(ttNIC_esa, details = TRUE, 
                       include = "?",
                       dir.exp = c(1, 1, 1, 1, 0, 0, 0))
# sol_Ni_esa

sol_Np_esa <- minimize(ttNIC_esa, include = "?", 
                       details = TRUE, show.cases = TRUE)

library(yarrr)
Ni_sol <- compute("RGOV*OPPG*~WEQU*~EUC + RGOV*~HRM*~WPOW*~WEQU*~EUC + LDB*RGOV*OPPG*HRM*~WPOW*~WEQU", data = ICdatay) # compute sol_Ni_esa set membership scores
col <- rep(yarrr::transparent("white", trans.val = .9), nrow(ICdatay)) # define a vector of colors
col[Ni_sol < 0.5 & ICdatay$IC < 0.5] <- "black" # zone 4
col[Ni_sol < 0.5 & ICdatay$IC >= 0.5] <- "grey60" # zone 5
col[Ni_sol > 0.5 & ICdatay$IC > 0.5 & ICdatay$IC < 0.9] <- "grey60" # zone 5
# col[nec > 0.5 & nec < 0.8 & ICdatay$IC > 0.5] <- "black"
clabels <- logical(nrow(ICdatay))
clabels[c(4, 7, 10, 18, 22, 31, 38, 2, 23, 33)] <- TRUE # only these labels to print

## Ni_sol
XYplot("RGOV*OPPG*~WEQU*~EUC + RGOV*~HRM*~WPOW*~WEQU*~EUC + LDB*RGOV*OPPG*HRM*~WPOW*~WEQU", ~IC, relation = "sufficiency", data = ICdatay, jitter = T, amount = 0.04, xlab = "Solution", clabels = clabels, cex = 0.7, col = "black", bg = col)


# Calibrated data ---------------------------------------------------------

CalData <- select(ICdatax, c("IC","LDB","RGOV","OPPG","HRM","WPOW","WEQU","EUC"))

rownames(CalData)[rownames(CalData) == "Bosnia and Herzegovina"] <- "BiH"
rownames(CalData)[rownames(CalData) == "North Macedonia"] <- "N. Macedonia"
rownames(CalData)[rownames(CalData) == "United Kingdom"] <- "U.K."

CalData %>% 
  kbl(caption = "Calibrated data. \\linebreak \\scriptsize{Key: OPPG = strong social opposition to gender equality, HRM = homogenous religious majority, WEQU = strong social support for women's equality, WPOW = strong women's empowerment, RGOV = solid right-wing government, LDB = liberal democratic backsliding, EUC = EU candidacy}. \\label{tab:calibrated-data}", digits=2) %>% 
  kable_classic_2(full_width = F, font_size = 9) %>% 
  row_spec(0, bold = T)


# Calibrated OPPG set and raw EVS and ESS data. ---------------------------

c1 <- select(ICdata, c("A124_09","D081_N","F118_N","hmsfmlsh","hmsacld","freehms"))
c1$OPPG <- ICdatax$OPPG
c1 <- c1[, c("OPPG","A124_09","D081_N","F118_N","hmsfmlsh","hmsacld","freehms")]

rownames(c1)[rownames(c1) == "Bosnia and Herzegovina"] <- "BiH"
rownames(c1)[rownames(c1) == "North Macedonia"] <- "N. Macedonia"
rownames(c1)[rownames(c1) == "United Kingdom"] <- "U.K."

c1 %>% 
  kbl(caption = "Calibrated OPPG set and raw EVS and ESS data. \\label{tab:c1}") %>% 
  kable_classic_2(full_width = F, font_size = 9) %>% 
  row_spec(0, bold = T)


# Calibrated HRM set and raw EVS data. ------------------------------------

c2 <- select(ICdata, c("F025_MAJ","F025_PERMAJ"))
c2$HRM <- ICdatax$HRM
c2 <- c2[, c("HRM","F025_MAJ","F025_PERMAJ")]

rownames(c2)[rownames(c2) == "Bosnia and Herzegovina"] <- "BiH"
rownames(c2)[rownames(c2) == "North Macedonia"] <- "N. Macedonia"
rownames(c2)[rownames(c2) == "United Kingdom"] <- "U.K."

# imputed values from Eurobarometer 2021: https://data.europa.eu/data/datasets/s2237_95_2_516_eng?locale=en
missing_Belgium <- c("Belgium")
c2 <- c2 %>% mutate(F025_MAJ = ifelse(row.names(c2) %in% missing_Belgium, 
                                      "Roman Catholic (44%)", F025_MAJ)) 
missing_Ireland <- c("Ireland")
c2 <- c2 %>% mutate(F025_MAJ = ifelse(row.names(c2) %in% missing_Ireland, 
                                      "Roman Catholic (54%)", F025_MAJ)) 
missing_Latvia <- c("Latvia")
c2 <- c2 %>% mutate(F025_MAJ = ifelse(row.names(c2) %in% missing_Latvia, 
                                      "Orthodox (28%)", F025_MAJ)) 

c2 %>% 
  kbl(caption = "Calibrated HRM set and raw EVS data. \\label{tab:c2}") %>%
  kable_classic_2(full_width = F, font_size = 9) %>% 
  row_spec(0, bold = T)

# Calibrated WEQU set and raw EVS and ESS data ----------------------------

c3 <- select(ICdata, c("C001","mnrgtjb"))
c3$WEQU <- ICdatax$WEQU
c3 <- c3[, c("WEQU","C001","mnrgtjb")]

rownames(c3)[rownames(c3) == "Bosnia and Herzegovina"] <- "BiH"
rownames(c3)[rownames(c3) == "North Macedonia"] <- "N. Macedonia"
rownames(c3)[rownames(c3) == "United Kingdom"] <- "U.K."

c3 %>% 
  kbl(caption = "Calibrated WEQU set and raw EVS and ESS data. \\label{tab:c3}") %>% 
  kable_classic_2(full_width = F, font_size = 9) %>% 
  row_spec(0, bold = T)

# Calibrated WPOW set and raw V-Dem data ----------------------------------

c4 <- select(ICdata, c("v2x_genpp_2020","v2x_cspart_2020"))
c4$WPOW <- ICdatax$WPOW
c4 <- c4[, c("WPOW","v2x_genpp_2020","v2x_cspart_2020")]

rownames(c4)[rownames(c4) == "Bosnia and Herzegovina"] <- "BiH"
rownames(c4)[rownames(c4) == "North Macedonia"] <- "N. Macedonia"
rownames(c4)[rownames(c4) == "United Kingdom"] <- "U.K."

c4 %>% 
  kbl(caption = "Calibrated WPOW set and raw V-Dem data. \\label{tab:c4}", digits = 2) %>% 
  kable_classic_2(full_width = F, font_size = 9) %>% 
  row_spec(0, bold = T)

# Democratic conditions ---------------------------------------------------

d1 <- select(ICdatax, c("LDB","v2x_libdem_2013","v2x_libdem_2021","EDB","SLD","SED"))

rownames(d1)[rownames(d1) == "Bosnia and Herzegovina"] <- "BiH"
rownames(d1)[rownames(d1) == "North Macedonia"] <- "N. Macedonia"
rownames(d1)[rownames(d1) == "United Kingdom"] <- "U.K."

d1 %>% 
  kbl(caption = "Democratic conditions: liberal democracy backsliding (LDB), raw V-Dem liberal democracy index scores from 2013 and 2021, as well as three alternative conditions, egalitarian democracy backsliding (EDB), strong liberal democracy (SLD), and strong egalitarian democracy (SED). \\label{tab:d1}", digits=2) %>% 
  kable_classic_2(full_width = F, font_size = 9) %>% 
  row_spec(0, bold = T)

# Conservative solution ---------------------------------------------------

stargazerSol(sol_c,
             outcome = "IC",
             show.cases = T,
             type = "latex",
             title = "",
             out = NULL,
             digits = 3)

# Intermediate solution ---------------------------------------------------

stargazerSol(sol_i,
             outcome = "IC",
             show.cases = T,
             type = "latex",
             title = "",
             out = NULL,
             digits = 3)

# Non-outcome analysis ----------------------------------------------------

SUIN_nIC <- superSubset(data = ICdatax,
                        conditions = c("LDB","RGOV","OPPG","HRM","WPOW","WEQU","EUC"),
                        outcome = "~IC",
                        relation = "nec",
                        incl.cut = 0.8,
                        cov.cut = 0.4,
                        ron.cut = 0.5,
                        depth = 2)
# SUIN_nIC

stargazerSol(sol_Nc,
             outcome = "~IC",
             show.cases = T,
             type = "latex",
             title = "",
             out = NULL,
             digits = 3)

stargazerSol(sol_Np,
             outcome = "~IC",
             show.cases = T,
             type = "latex",
             title = "",
             out = NULL,
             digits = 3)

# QCA robustness protocol -------------------------------------------------

rob.calibrange(raw.data = ICdatax,
               calib.data = ICdatay,
               test.cond.raw = "LDEMBACK",
               test.cond.calib = "LDB",
               test.thresholds = c(0, -0.08, -0.1),
               type = "fuzzy",
               step = 0.01,
               max.runs = 20,
               outcome = "IC",
               conditions = c("LDB","RGOV","OPPG","HRM","WPOW","WEQU","EUC"),
               incl.cut = 0.85,
               include = "?")

rob.calibrange(raw.data = ICdatax,
               calib.data = ICdatay,
               test.cond.raw = "WOMPOW",
               test.cond.calib = "WPOW",
               test.thresholds = c(0.35, 0.9, 0.97),
               type = "fuzzy",
               step = 0.01,
               max.runs = 20,
               outcome = "IC",
               conditions = c("LDB","RGOV","OPPG","HRM","WPOW","WEQU","EUC"),
               incl.cut=0.85,
               include = "?")

# consistency and frequency -----------------------------------------

rob.inclrange(data = ICdatay,
              step = 0.01,
              max.runs = 20,
              outcome = "IC",
              conditions = c("LDB","RGOV","OPPG","HRM","WPOW","WEQU","EUC"),
              incl.cut=0.85,
              include = "?")

rob.ncutrange(data = ICdatay,
              step = 1,
              max.runs = 10,
              outcome = "IC",
              conditions = c("LDB","RGOV","OPPG","HRM","WPOW","WEQU","EUC"),
              incl.cut=0.85,
              n.cut = 1,
              include = "?")

# TS1 (changes to include inconsistent row in minimisation)

ttICr1 <- truthTable(data=ICdatay, outcome = "IC",
                     conditions = c("LDB","RGOV","OPPG","HRM","WPOW","WEQU","EUC"),
                     incl.cut=0.8,
                     sort.by=c("incl", "n"), 
                     complete=TRUE,
                     show.cases=TRUE, 
                     PRI=TRUE)
ttICr1

sol_pr1 <- minimize(ttICr1, include = "?", 
                    details = TRUE, show.cases = TRUE)
sol_pr1

# TS2a 

# LDB: liberal democratic backsliding
ICdatax$LDB <- calibrate(ICdatax$LDEMBACK,
                         type = "fuzzy",
                         thresholds = c(0.02, -0.11, -0.25))

ICdataR2A <- select(ICdatax, c("IC","LDB","RGOV","OPPG","HRM","WPOW","WEQU","EUC"))

ttICr2a <- truthTable(data=ICdataR2A, outcome = "IC",
                      conditions = c("LDB","RGOV","OPPG","HRM","WPOW","WEQU","EUC"),
                      incl.cut=0.8,
                      sort.by=c("incl", "n"), 
                      complete=TRUE,
                      show.cases=TRUE, 
                      PRI=TRUE)
ttICr2a

sol_pr2a <- minimize(ttICr2a, include = "?", 
                     details = TRUE, show.cases = TRUE)
sol_pr2a

# restore original calibration of LDB
ICdatax$LDB <- calibrate(ICdatax$LDEMBACK,
                         type = "fuzzy",
                         thresholds = c(0, -0.08, -0.1))

# TS2b

# WPOW: women's political empowerment 
ICdatax$WPOW <- calibrate(ICdatax$WOMPOW,
                          type = "fuzzy",
                          thresholds = c(0.3, 0.85, 0.93))

ICdataR2B <- select(ICdatax, c("IC","LDB","RGOV","OPPG","HRM","WPOW","WEQU","EUC"))

ttICr2b <- truthTable(data=ICdataR2B, outcome = "IC",
                      conditions = c("LDB","RGOV","OPPG","HRM","WPOW","WEQU","EUC"),
                      incl.cut=0.8,
                      sort.by=c("incl", "n"), 
                      complete=TRUE,
                      show.cases=TRUE, 
                      PRI=TRUE)
ttICr2b

sol_pr2b <- minimize(ttICr2b, include = "?", 
                     details = TRUE, show.cases = TRUE)
sol_pr2b

# Obtaining test solution (TS) and robust core (RC) 

TS <- list(sol_pr1,sol_pr2a,sol_pr2b)

# Obtain parameters of fit for RC:
rob.corefit(test_sol = TS,
            initial_sol = sol_p,
            outcome = "IC")

# Obtain Boolean expression for RC:
intersection(sol_p, sol_pr1, sol_pr2a, sol_pr2b)

# Identify Robustness-relevant Types of Cases and the RCRs
o <- rob.xyplot(test_sol = TS,
                initial_sol = sol_p,
                outcome = "IC",
                all_labels = TRUE,
                fontsize = 3,
                jitter = TRUE,
                area_lab = TRUE)
o

rob.cases(test_sol = TS,
          initial_sol = sol_p,
          outcome = "IC")

rob.singletest(test_sol = TS,
               initial_sol = sol_p,
               outcome = "IC")

# Alternative analyses, different democracy measures -----------

### Necessity analysis (EDB)

ICdataEDB <- select(ICdatax, c("IC","EDB","RGOV","OPPG","HRM","WPOW","WEQU","EUC"))

SUIN_ICedb <- superSubset(data = ICdataEDB,
                          conditions = c("EDB","RGOV","OPPG","HRM","WPOW","WEQU","EUC"),
                          outcome = "IC",
                          relation = "nec",
                          incl.cut = 0.9,
                          cov.cut = 0.6,
                          ron.cut = 0.5,
                          depth = 2)
# SUIN_ICedb

# pimplot(data = ICdataEDB, results = SUIN_ICedb, outcome = "IC",
#         necessity = TRUE, jitter = TRUE, font = "mono")

### Sufficiency analysis (EDB)

# create a truth table
ttIC2 <- truthTable(data=ICdataEDB, outcome = "IC",
                    conditions = c("EDB","RGOV","OPPG","HRM","WPOW","WEQU","EUC"),
                    incl.cut=0.9,
                    sort.by=c("incl", "n"), 
                    complete=TRUE,
                    show.cases=TRUE, 
                    PRI=TRUE)
# ttIC2

# pimplot(data=ttIC2, results=sol_c, ttrows=c(79), outcome="IC", all_labels = T, jitter=T)

sol_c2 <- minimize(ttIC2, details=TRUE, show.cases=TRUE)
# sol_c2
# model ambiguity

# stargazerSol(sol_c2,
#              outcome = "IC",
#              sol = 1,
#              show.cases = T,
#              type = "latex", 
#              title = "",
#              out = NULL,
#              digits = 3)

# pimplot(data = ttIC2, results = sol_c2, outcome = "IC", all_labels = TRUE, jitter = TRUE, font = "mono")


sol_i2 <- minimize(ttIC2, details = TRUE, 
                   include = "?",
                   dir.exp = c(0, 0, 0, 0, 1, 1, 1))
# sol_i2

# stargazerSol(sol_i2,
#              outcome = "IC",
#              sol = 1,
#              show.cases = T,
#              type = "latex",
#              title = "",
#              out = NULL,
#              digits = 3)

# pimplot(data = ttIC2, results = sol_i2, outcome = "IC",
#         all_labels = TRUE, jitter = TRUE, font = "mono")

sol_p2 <- minimize(ttIC2, include = "?", 
                   details = TRUE, show.cases = TRUE)
# sol_p2
# model ambiguity

# stargazerSol(sol_p2,
#              outcome = "IC",
#              sol = 1,
#              show.cases = T,
#              type = "latex", 
#              title = "",
#              out = NULL,
#              digits = 3)

# pimplot(data = ttIC2, results = sol_p2, outcome = "IC", all_labels = TRUE, jitter = TRUE, font = "mono")

# smmr(results = sol_p2, outcome = "IC", match=FALSE, cases=6, term = 2)

## Strong liberal democracy (SLD)

### Necessity analysis (SLD)

ICdataSLD <- select(ICdatax, c("IC","SLD","RGOV","OPPG","HRM","WPOW","WEQU","EUC"))

SUIN_ICSLD <- superSubset(data = ICdataSLD,
                          conditions = c("SLD","RGOV","OPPG","HRM","WPOW","WEQU","EUC"),
                          outcome = "IC",
                          relation = "nec",
                          incl.cut = 0.9,
                          cov.cut = 0.6,
                          ron.cut = 0.5,
                          depth = 2)
# SUIN_ICSLD

# pimplot(data = ICdataSLD, results = SUIN_ICSLD, outcome = "IC",
#         necessity = TRUE, jitter = TRUE, font = "mono")

### Sufficiency analysis (SLD)

# create a truth table
ttIC2 <- truthTable(data=ICdataSLD, outcome = "IC",
                    conditions = c("SLD","RGOV","OPPG","HRM","WPOW","WEQU","EUC"),
                    incl.cut=0.9,
                    sort.by=c("incl", "n"), 
                    complete=TRUE,
                    show.cases=TRUE, 
                    PRI=TRUE)
# ttIC2

# pimplot(data=ttIC2, results=sol_c, ttrows=c(79), outcome="IC", all_labels = T, jitter=T)

sol_c2 <- minimize(ttIC2, details=TRUE, show.cases=TRUE)
# sol_c2
# model ambiguity

# stargazerSol(sol_c2,
#              outcome = "IC",
#              sol = 1,
#              show.cases = T,
#              type = "latex",
#              title = "",
#              out = NULL,
#              digits = 3)

# pimplot(data = ttIC2, results = sol_c2, outcome = "IC", all_labels = TRUE, jitter = TRUE, font = "mono")

sol_i2 <- minimize(ttIC2, details = TRUE, 
                   include = "?",
                   dir.exp = c(0, 0, 0, 0, 1, 1, 1))
# sol_i2

# stargazerSol(sol_i2,
#              outcome = "IC",
#              sol = 1,
#              show.cases = T,
#              type = "latex",
#              title = "",
#              out = NULL,
#              digits = 3)

# pimplot(data = ttIC2, results = sol_i2, outcome = "IC",
#         all_labels = TRUE, jitter = TRUE, font = "mono")


sol_p2 <- minimize(ttIC2, include = "?", 
                   details = TRUE, show.cases = TRUE)
# sol_p2
# model ambiguity

# stargazerSol(sol_p2,
#              outcome = "IC",
#              sol = 1,
#              show.cases = T,
#              type = "latex", 
#              title = "",
#              out = NULL,
#              digits = 3)

# pimplot(data = ttIC2, results = sol_p2, outcome = "IC", all_labels = TRUE, jitter = TRUE, font = "mono")

# smmr(results = sol_p2, outcome = "IC", match=FALSE, cases=6, term = 2)

## Strong egalitarian democracy (SED)

### Necessity analysis (SED)

ICdataSED <- select(ICdatax, c("IC","SED","RGOV","OPPG","HRM","WPOW","WEQU","EUC"))

SUIN_ICSED <- superSubset(data = ICdataSED,
                          conditions = c("SED","RGOV","OPPG","HRM","WPOW","WEQU","EUC"),
                          outcome = "IC",
                          relation = "nec",
                          incl.cut = 0.9,
                          cov.cut = 0.6,
                          ron.cut = 0.5,
                          depth = 2)
# SUIN_ICSED

# pimplot(data = ICdataSED, results = SUIN_ICSED, outcome = "IC",
#         necessity = TRUE, jitter = TRUE, font = "mono")

### Sufficiency analysis (SED)

# create a truth table
ttIC2 <- truthTable(data=ICdataSED, outcome = "IC",
                    conditions = c("SED","RGOV","OPPG","HRM","WPOW","WEQU","EUC"),
                    incl.cut=0.9,
                    sort.by=c("incl", "n"), 
                    complete=TRUE,
                    show.cases=TRUE, 
                    PRI=TRUE)
# ttIC2

# pimplot(data=ttIC2, results=sol_c, ttrows=c(79), outcome="IC", all_labels = T, jitter=T)

sol_c2 <- minimize(ttIC2, details=TRUE, show.cases=TRUE)
# sol_c2
# model ambiguity

# stargazerSol(sol_c2,
#              outcome = "IC",
#              sol = 1,
#              show.cases = T,
#              type = "latex",
#              title = "",
#              out = NULL,
#              digits = 3)

# pimplot(data = ttIC2, results = sol_c2, outcome = "IC", all_labels = TRUE, jitter = TRUE, font = "mono")

sol_i2 <- minimize(ttIC2, details = TRUE, 
                   include = "?",
                   dir.exp = c(0, 0, 0, 0, 1, 1, 1))
# sol_i2

# stargazerSol(sol_i2,
#              outcome = "IC",
#              sol = 1,
#              show.cases = T,
#              type = "latex",
#              title = "",
#              out = NULL,
#              digits = 3)

# pimplot(data = ttIC2, results = sol_i2, outcome = "IC",
#         all_labels = TRUE, jitter = TRUE, font = "mono")

sol_p2 <- minimize(ttIC2, include = "?", 
                   details = TRUE, show.cases = TRUE)
# sol_p2
# model ambiguity

# stargazerSol(sol_p2,
#              outcome = "IC",
#              sol = 1,
#              show.cases = T,
#              type = "latex",
#              title = "",
#              out = NULL,
#              digits = 3)

# pimplot(data = ttIC2, results = sol_p2, outcome = "IC", all_labels = TRUE, jitter = TRUE, font = "mono")

# smmr(results = sol_p2, outcome = "IC", match=FALSE, cases=6, term = 2)

